/**
 * @file Multigrid.h
 * @author Lishijie (lsj1018845759@outlook.com)
 * @brief 
 * @version 0.1
 * @date 2020-10-14
 * 
 * @copyright Copyright (c) 2020
 * 
 */
#ifndef _Multigrid_H_
#define _Multigrid_H_

#include "Restrict.h"
#include "Interpolation.h"
#include <vector>
#include <algorithm>
#include <stdio.h>
#include <iostream>
#include <string>
#include <cstring>
#include <cmath>
#include <math.h>
#include "Miscellaneous.h"
double norm_max(std::vector<double>);
double norm_2(std::vector<double>);
namespace MyMultigrid
{
    /**
     * @brief The options of this Multigrid solver
     *       PROVIDEDOPTIONS:  
     *       1:BoundoryConditionType:Homogeneous(default) OR Nohomogeneous;
     *       2:RestrictionOperatorType:Full-weighting(default) OR Injection;
     *       3:InterpolationOperatorOperatorType:Linear(default) OR Quadratic;
     *       4:CycleType:Vcycle(default) OR FMcycle;
     *       5:MaxIterNum:Integer(default 10);
     *       6:RelAccurcy:double(default 1e - 6);
     */ 
    struct Options
    {
        std::string BoundoryConditionType = "Homogeneous";
        std::string RestrictionOperatorType = "Full-weighting";
        std::string InterpolationOperatorOperatorType = "Linear";
        std::string CycleType = "Vcycle";
        int MaxIterNum = 10;
        double RelAccurcy = 1e-5;
    };

    #define TEMPLATE template<int N>
    TEMPLATE
    /**
     * @brief Template Multigrid solver
     *        N: The Number of segments in interval [0 , 1];
     */
    class Multigrid
    {
    private:
        
        /**
         * @brief Storage the variable of v or u or e;
         */
        std::vector<double> v;
        /**
         * @brief Storage the variable of f or r;
         */
        std::vector<double> f;
        /**
         * @brief Storage the solution solved by Mutigrid method for each iteration;
         */
        std::vector<double> solution;
        /**
         * @brief Storage the  true solution of this equation; 
         */
        std::vector<double> real_solution;
        /**
         * @brief Set the option of this Mutigrid method;
         */
        Options option;
        /**
         * @brief Storage the real solution function of this problem;
         *         Pointer to function;
         */
        double(*real)(double);
        /**
         * @brief Storage the righthands terms function of possion equation
         *        Pointer to function;
         */
        double(*func)(double);
        /**
         * @brief The highest level of the N;
         */
        
        int level =  int(log(N) / log(2));
        /**
         * @brief Selection of restrict operator
         */
        RestrictOperator *restrict;
        /**
         * @brief Selection of Interpolation operator
         */
        InterpolationOperator *interpolation;
        /**
         * @brief current relativeaccuracy;
         */
        double relativeaccuracy;
        /**
         * @brief current iteration number
         */
        int IterNumber = 0;
        /**
         * @brief Dirichlet boundary condition
         */
        double u0 , u1;

    public:
        /**
         * @brief Construct a new Multigrid object
         */
        Multigrid()
        {}; 
        /**
         * @brief Construct a new Multigrid object
         * 
         * @param u real soution function of possion equation
         * @param fun righthands term function of possion equation
         * @param opt The options of Multigrid method
         */
        Multigrid(double(*u)(double),double(*fun)(double),Options opt);
        /**
         * @brief Destroy the Multigrid object
         */
        ~Multigrid(){};
        /**
         * @brief solve this problem by Multigrid method
         */
        void solve();
        /**
         * @brief overload-output operater
         * 
         * @tparam N_ 
         * @param os 
         * @param m 
         * @return std::ostream& 
         */
        template<int N_>
        friend std::ostream& operator<<(std::ostream & os,Multigrid<N_> m);

    private:
        /**
         * @brief initialize or reset the Multigrid method;
         */
        void init();

        /**
         * @brief The Weighted-jacobi iteration method
         * 
         * @param l current level
         * @param times iteration times
         * @param omage weight
         */
        void weightedjacobi(int l,int times,double omage,int s);
        
        
        /**
         * @brief Vcycle iteration strategy
         * 
         * @param l current level
         */
        void Vcycle(int l,int s) ;
        /**
         * @brief FMcycle iteration strategy
         * 
         * @param l current level
         */
        void FMcycle(int l);
        /**
         * @brief Construct a new operator<< object
         * 
         * @tparam N_ he Number of segments in interval [0 , 1];
         * @param os out ostream
         * @param m  Multigrid solver
         */
        int low(int l){return pow(2,l) + l - 2;};
        int high(int l){return pow(2,l + 1) + l - 2;};
        void resetv(int l){ for(int i = 0;i <= high(l);i++) v[i] = 0;};   
    };

    TEMPLATE
    Multigrid<N>::Multigrid(double(*u)(double),double(*fun)(double) ,Options opt)
    {
        
        v.resize(pow(2,level + 1) + level - 1);
        f.resize(pow(2,level + 1) + level - 1);
        int lo = low(level);
        int hi = high(level);
        double h = 1 / double(pow(2,level));     
        for(int i = lo ;i <= hi;i++)
        {
            f[i] = fun(h * (i - lo));
            real_solution.push_back(u(h * (i - lo)));
        }
        solution.resize(real_solution.size());
        real = u;
        func = fun;
        option = opt;
        u0 = real(0);
        u1 = real(1);
        if(opt.RestrictionOperatorType == "Full-weighting")
            restrict = new FullWeightingRestrictOperator(v);
        else
        {
            restrict = new InjectionRestrictOperator(v);
        }

        if(opt.InterpolationOperatorOperatorType == "Linear")
            interpolation = new LinearInterpolaionOperator(v);
        else
        {
            interpolation = new QuadraticInterpolaionOperator(v);
        }
        
    };

    TEMPLATE
    void Multigrid<N>::init()
    {
        f.clear();
        f.resize(pow(2,level + 1) + level - 1);
        int lo = low(level);
        int hi = high(level);
        f[lo] = u0;
        f[hi] = u1;
        double h = 1 / double(pow(2,level));
        for(int i = lo + 1;i < hi;i++)
            f[i] = func(h * (i - lo));
        for(int i = 1;i < lo;i++)
            v[i] = 0.0;
    };

    TEMPLATE
    void Multigrid<N>::Vcycle(int l,int s)
    {
        std::vector<double> r;
        r = f;
        weightedjacobi(l,10,2/3.0,s);
        if(l > 1)
        {
            int lo = low(l);
            int hi = high(l); 
            r[lo] = 0;
            r[hi] = 0;
            for(int i = low(l) + 1;i < high(l);i++)
            { 
                double h = 1.0 /(double(pow(2,l)));
                r[i] = f[i] - ((-v[i - 1] + 2 * v[i] -v[i + 1])/(h * h));
            }
            restrict -> set(r);
            restrict -> do_restrict(l);
            std::copy(restrict -> input.begin() + low(l - 1),restrict -> input.begin() + high(l -1),f.begin()+ low(l - 1));
            Vcycle(l-1,s);
        }
        else
        {
            int current = low(l) + 1;
            double h = 1.0 /(double(pow(2,l)));
            v[current] = 0.5*(h * h * f[current]);//高斯消去直接精确求解
        }
        if( l != s)
        {
            std::vector<double> t = v;
            t[low(l)] = 0;
            t[high(l)] = 0;
            interpolation->set(t);
            interpolation->do_Interpolation(l);
            for(int i = low(l + 1) + 1;i < high(l + 1);i++)
            {
                v[i] = interpolation->input[i] + v[i];
            }
            weightedjacobi(l + 1,1,2/3.0,s);
        }

    };
    TEMPLATE
    void Multigrid<N>::FMcycle(int l)
    {
        if(l > 1)
        {
            
            int lo = low(l);
            int hi = high(l);
            restrict->set(f);
            restrict->do_restrict(l);
            std::copy(restrict -> input.begin() + low(l - 1),restrict -> input.begin() + high(l -1) + 1,f.begin()+ low(l - 1));
            FMcycle(l-1);
        }
        else
        {
            int current = low(l) + 1;
            double h = 1.0 /(double(pow(2,l)));
            v[current] = 0.5*(h * h * f[current] + u0 + u1);
        }
        if(l != level)
        {
            interpolation->set(v);
            interpolation->do_Interpolation(l);
            v.erase(v.begin() + low(l + 1),v.begin() +high(l + 1));
            v.insert(v.begin() + low(l + 1),interpolation ->input.begin()+low(l + 1),interpolation -> input.begin()+high(l + 1));       
            for(int i = 0;i < 2;i++)
            {
                resetv(l);
                Vcycle(l + 1,l + 1);
            }    
        }  
    }
    
    TEMPLATE
    void Multigrid<N>::solve()
    {
        std::vector<double> temp(real_solution.size());
        for(int i = 0;i < temp.size();i++)
        {
            temp[i] = real_solution[i] - solution[i];
        }
        double acctemp = norm_max(temp);
        int currenttimes = 0;
        int lo = low(level);
        int hi = high(level);
        std::cout << "..Calculating.."<<std::endl;
        do
        {
            
            if(option.CycleType == "Vcycle")
                Vcycle(level,level);
            else
                FMcycle(level);
            for(int i = lo;i <= hi;i++)
            {
                solution[i - lo] = v[i];
            }
            for(int i = 0;i < temp.size();i++)
            {
                temp[i] = real_solution[i] - solution[i];
            }
            acctemp = norm_max(temp);
            currenttimes++;
            init();
            std::cout <<"Current iteration time is :" << currenttimes << std::endl;
            printf("Current relative accurcy is : %.8lf\n",acctemp);
        }while (currenttimes < option.MaxIterNum && acctemp > option.RelAccurcy);
        relativeaccuracy = acctemp;
        IterNumber = currenttimes;
        if(currenttimes >= option.MaxIterNum && relativeaccuracy > option.RelAccurcy)
        {
            std::cout << "Iteration number arrive the Max Number!"<<std::endl;
            std::cout << "But does not convergence to the preset relative accrucy!"<<std::endl;
        }
        else
        {
            std::cout << "Iteration convergence to " << option.RelAccurcy << " at "<< IterNumber << "th time"<<std::endl;
            printf("The last relative accurcy is : %.8lf !",relativeaccuracy);
        } 
        std::cout << "..Done!.."<<std::endl;
    }; 

    TEMPLATE
    void Multigrid<N>::weightedjacobi(int l,int times,double omage,int s)
    {
        std::vector<double> temp;
        if(l == s)
        {
            v[low(l)] = u0;
            v[high(l)] = u1;
        }
        else
        {
            v[low(l)] = 0;
            v[high(l)] = 0;
        }
        for(unsigned int k = 0;k < times;k++)
        {
            temp = v;
            for(int i = low(l) + 1;i <= high(l) - 1;i++)
            {
                double t = 0.5*(temp[i - 1] + temp[i + 1] + (1/double(pow(2,l))) * (1/double(pow(2,l)))*f[i]);
                v[i] = temp[i] + omage *(t - temp[i]);
            }
        }
    }
    template<int N_>
    std::ostream& operator<<(std::ostream& os,Multigrid<N_> m)
    {
        os << "MyMultigrid::Multigrid Options setting information !" <<std::endl;
        os << "The segments of the Multigrid solver : " << N_ <<std::endl;
        os << "The Options of the mutigrid solver : " << m.option.BoundoryConditionType <<std::endl;
        os << "                                     " <<m.option.RestrictionOperatorType << std::endl;
        os << "                                     " << m.option.InterpolationOperatorOperatorType << std::endl;
        os << "                                     " << m.option.CycleType<<std::endl;
        os << "                                     " << m.option.MaxIterNum<<std::endl;
        os << "                                     " <<m.option.RelAccurcy <<std::endl;
        return os;
    }
} // namespace MyMultigrid

    

#else

#endif


